home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Best of www.BestZips.com (Collector's Edition)
/
Best of WWW.BESTZIPS.COM Collector's Edition (JCSM Shareware) (JCS Marketing).ISO
/
prgtools
/
tn2.zip
/
SUNUP.T
< prev
next >
Wrap
Text File
|
1996-11-15
|
10KB
|
442 lines
%
% "sunup.t" calculates the times of sunrise and sunset
% on any date, accurate to the minute within several
% centuries of the present. It correctly describes what
% happens in the arctic and antarctic regions, where
% the Sun may not rise or set on a given date.
%
% This program was adapted from a BASIC program in
% Sky & Telescope magazine, August 1994, page 84,
% written by Roger W. Sinnott.
%
% Enter north latitudes positive, west longitudes negative.
%
% For the time zone, enter the number of hours west of Greenwich
% (e.g., 5 for EST, 4 for EDT).
%
% Sample program for the T Interpreter by:
%
% Stephen R. Schmitt
% 962 Depot Road
% Boxborough, MA 01719
%
const PI : real := 2 * arcsin( 1 )
const DR : real := PI / 180
const K1 : real := 15 * DR * 1.0027379
var Sunrise, Sunset : boolean := false
type rvector : array[3] of real
program
var lat, lon, tz, ph : real
var jd, ct, t0, ra0, ra1, dec0, dec1 : real
var k, zone : int
var ra, dec, v : rvector
prompt "Latitude, +N, -S (deg):"
get lat
prompt "Longitude, +E, -W (deg):"
get lon
prompt "Time zone +West of GMT, -East (hrs):"
get zone
jd := julian_day - 2451545 % Julian day relative to Jan 1.5, 2000
show_input( lat, lon, zone )
lon := lon / 360
tz := zone / 24
ct := jd / 36525 + 1 % centuries since 1900.0
t0 := lst( lon, jd, tz )
% get sun's position
jd := jd + tz
sun( ra0, dec0, jd, ct )
jd := jd + 1
sun( ra1, dec1, jd, ct )
if ra1 < ra0 then
ra1 := ra1 + 2 * PI
end if
ra[0] := ra0
dec[0] := dec0
for k := 0...23 do
ph := ( k + 1 ) / 24
ra[2] := ra0 + ph * ( ra1 - ra0 )
dec[2] := dec0 + ph * ( dec1 - dec0 )
test( k, zone, t0, lat, ra, dec, v )
ra[0] := ra[2]
dec[0] := dec[2]
end for
% special message?
if ( not Sunrise ) and ( not Sunset ) then
if v[2] < 0 then
put "Sun down all day"
else
put "Sun up all day"
end if
else
if not Sunrise then
put "No sunrise this date"
elsif not Sunset then
put "No sunset this date"
end if
end if
end program
%
% display latitude, longitude and time zone
%
procedure show_input( lat, lon : real, zone : int )
var deg, min : int
if sgn( zone ) = sgn( lon ) and zone ~= 0 then
put "WARNING: time zone and longitude are incompatible"
end if
if abs( zone ) > 12 then
put "WARNING: invalid time zone"
end if
if abs( lat ) > 90 then
put "WARNING: invalid latitude"
end if
if abs( lon ) > 180 then
put "WARNING: invalid longitude"
end if
if lat > 0.0 then
deg := floor( lat )
min := floor( 60 * ( lat - deg ) )
put deg:4, ":",min:2, " N "...
else
deg := floor( -lat )
min := floor( 60 * ( -lat - deg ) )
put deg:4, ":",min:2, " S "...
end if
if lon > 0.0 then
deg := floor( lon )
min := floor( 60 * ( lon - deg ) )
put deg:4, ":",min:2, " E"
else
deg := floor( -lon )
min := floor( 60 * ( -lon - deg ) )
put deg:4, ":",min:2, " W"
end if
end procedure
%
% LST at 0h zone time
%
function lst( lon, jd, z : real ) : real
var s, t0 : real
s := 24110.5 + 8640184.812999999 * jd / 36525
s := s + 86636.6 * z + 86400 * lon
s := s / 86400
s := s - floor( s )
t0 := s * 360 * DR
return t0
end function
%
% test an hour for an event
%
procedure test( k, zone : int,
t0, lat : real,
var ra, dec, v : rvector )
var ha : rvector
var a, b, c, d, e, s, z : real
var hr, min, time : real
var az, dz, hz, nz : real
var zchar : string
var zlet : string := "MLKIHGFEDCBAZNOPQRSTUVWXY"
label test_exit :
ha[0] := t0 - ra[0] + k * K1
ha[2] := t0 - ra[2] + k * K1 + K1
ha[1] := ( ha[2] + ha[0] ) / 2 % hour angle at half hour
dec[1] := ( dec[2] + dec[0] ) / 2 % declination at half hour
s := sin( lat * DR )
c := cos( lat * DR )
z := cos( 90.833 * DR )
if k <= 0 then
v[0] := s * sin( dec[0] ) + c * cos( dec[0] ) * cos( ha[0] ) - z
end if
v[2] := s * sin( dec[2] ) + c * cos( dec[2] ) * cos( ha[2] ) - z
if sgn( v[0] ) = sgn( v[2] ) then
goto test_exit
end if
v[1] := s * sin( dec[1] ) + c * cos( dec[1] ) * cos( ha[1] ) - z
a := 2 * v[0] - 4 * v[1] + 2 * v[2]
b := -3 * v[0] + 4 * v[1] - v[2]
d := b * b - 4 * a * v[0]
if d < 0 then
goto test_exit
end if
d := sqrt( d )
if v[0] < 0 and v[2] > 0 then
put "Sunrise at: "...
Sunrise := true
end if
if v[0] > 0 and v[2] < 0 then
put "Sunset at: "...
Sunset := true
end if
e := ( -b + d ) / ( 2 * a )
if e > 1 or e < 0 then
e := ( -b - d ) / ( 2 * a )
end if
% time of an event
time := k + e + 1 / 120
hr := floor( time )
min := floor( ( time - hr ) * 60 )
zchar := " "
if abs( zone ) <= 12 then
zchar[1] := zlet[zone+12]
end if
put hr:2, ":", min:2, zchar...
% azimuth of the sun at the event
hz := ha[0] + e * ( ha[2] - ha[0] )
nz := -cos( dec[1] ) * sin( hz )
dz := c * sin( dec[1] ) - s * cos( dec[1] ) * cos( hz )
az := arctan( nz / dz ) / DR
if dz < 0 then
az := az + 180
end if
if az < 0 then
az := az + 360
end if
if az > 360 then
az := az - 360
end if
put "azimuth: ", az:5:1
test_exit:
v[0] := v[2]
end procedure
%
% sun's position using fundamental arguments
% (Van Flandern & Pulkkinen, 1979)
%
procedure sun( var ra, dec : real, jd, ct : real )
var g, lo, s, u, v, w : real
lo := 0.779072 + 0.00273790931 * jd
lo := lo - floor( lo )
lo := lo * 2 * PI
g := 0.993126 + 0.0027377785 * jd
g := g - floor( g )
g := g * 2 * PI
v := 0.39785 * sin( lo )
v := v - 0.01 * sin( lo - g )
v := v + 0.00333 * sin( lo + g )
v := v - 0.00021 * ct * sin( lo )
u := 1 - 0.03349 * cos( g )
u := u - 0.00014 * cos( 2 * lo )
u := u + 0.00008 * cos( lo )
w := -0.0001 - 0.04129 * sin( 2 * lo )
w := w + 0.03211 * sin( g )
w := w + 0.00104 * sin( 2 * lo - g )
w := w - 0.00035 * sin( 2 * lo + g )
w := w - 0.00008 * ct * sin( g )
% compute sun's right ascension and declination
s := w / sqrt( u - v * v )
ra := lo + arctan( s / sqrt( 1 - s * s ) )
s := v / sqrt( u )
dec := arctan( s / sqrt( 1 - s * s ) )
end procedure
%
% determine Julian day from calendar date
% ref: Jean Meeus, "Astronomical Algorithms", Willmann-Bell, 1991
%
function julian_day : real
var a, b, day, jd : real
var month, year : int
var gregorian : boolean
prompt "Year:"
get year
prompt "Month:"
get month
prompt "Day:"
get day
date( floor( day ), month, year )
if year < 1583 then
gregorian := false
else
gregorian := true
end if
if month = 1 or month = 2 then
year := year - 1
month := month + 12
end if
a := floor( year / 100 )
if gregorian then
b := 2 - a + floor( a / 4 )
else
b := 0.0
end if
jd := floor( 365.25 * ( year + 4716 ) ) +
floor( 30.6001 * ( month + 1 ) ) +
day + b - 1524.5
return jd
end function
%
% display the date
%
procedure date( d, m, y : int )
var day, month : array[12] of string
var mi, yi, i : int
month[0] := "January"
month[1] := "February"
month[2] := "March"
month[3] := "April"
month[4] := "May"
month[5] := "June"
month[6] := "July"
month[7] := "August"
month[8] := "September"
month[9] := "October"
month[10] := "November"
month[11] := "December"
day[0] := "Monday"
day[1] := "Tuesday"
day[2] := "Wednesday"
day[3] := "Thursday"
day[4] := "Friday"
day[5] := "Saturday"
day[6] := "Sunday"
mi := m
yi := y
if y < 1752 then
put month[mi-1], " ", d, ", ", yi
else
if m = 1 or m = 2 then
m := m + 12
y := y - 1
end if
i := d + 2*m + 3*(m+1) div 5 + y + y div 4 - y div 100 + y div 400
i := i mod 7
put day[i], ", ", month[mi-1], " ", d, ", ", yi
end if
end procedure
%
% returns value for sign of argument
%
function sgn( x : real ) : int
var rv : int
if x > 0.0 then
rv := 1
elsif x < 0.0 then
rv := -1
else
rv := 0
end if
return rv
end function
%
% absolute value of argument
%
function abs( x : real ) : real
if x < 0.0 then
x := -x
end if
return x
end function